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ABSTRACT 


A finite element formulation for solving axisymmetric 
meansitent heat conduction and thermal stress problems is 
developed in this thesis. The governing equations of 
EDoUu»ted, lanear, isotropic thermoelasticity are discretized 
Hsing Quadratic isoparametric elements. A FORTRAN IV program, 
using double precision arithmetic, is presented. Compact 
storage techniques for banded symmetric matrices are used. 

Comparisons between exact and computer solutions demon- 
EM е Е1оѕе agreement for a number of test problems. De- 


talled instructions for using the program are included. 
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Note: 


ti | o Ia lo Io  οὁ 


O 


E τ 15 [ж de 1н 


LISI UOT ΝΟ 


ΙΙ; πε 1: used to denote a column vector 
and a double underline denotes a rectangular matrix. 
The symbols used in computer program are described in 
the beginning of Appendix B. 

Entrance fluid cross-sectional area 

Base recombination of C and Y matrices 

A constant vector 

Standard rectangular strain-displacement matrix 
Thermal capacitance matrix 

Specific heat 

Standard elasticity matrix 

Young's modulus of elasticity 
ENDcrscrIdMDIdesTonatamp element contribution 

Load vector 

Element thermal load vector 

Element pressure load vector 

Element centrifugal load vector 

neat combination Ol € and Y matrices 

Surface heat transfer coefficient 

Identity matrix 

Jacobian coordinate transformation matrix 

System stiffness matrix 

Element stiffness matrix 

Thermal conductivity 


Thickness of slab 





с [C - ы - [ч = wa әу з 


lz = < I< 


со чс A — < xX K Z ες 
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Po 
lo 
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=) 
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аур 
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Arc length along the side of quadrilateral 
Shape function 

Total number of nodes 

Outward normal or number of nodes per element 
Pressure 

Radial coordinate 

Surface area 

Nodal temperature vector 


Temperature or, when used as a superscript, 
transpose of a matrix 


Average temperature 

Fluid temperature (used in one-dimensional example) 
A vector defined as <1 1 1 o>! 
Radial displacement 
Right-hand side vector in conduction equation 
Volume 

Work done by loads 

Modal matrix 

Axial displacement 

Eigenvector 

Thermal admittance matrix 

Dependent variable 

ΠΕ ij Of the matrix Y* 

Matrix representation 

Row vector 

Gradient operator 

Coefficient of thermal expansion 


Constant coefficient vector 
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Nodal displacement vector 
Element displacement vector 
Straln vector 


Eigenvalue of one-dimensional transient temperature 
solution 


Thermal strain vector 
Element strain vector 

Local element coordinate 
Poisson's ratio 

Eigenvalue 

Material density 

Stress vector 

Local element coordinate 
Time 

Step size of numerical time integration 
Shearing stress component 
Fluid temperature, or angle 
Spectral matrix 


Speed of rotation 
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TS INIRODUCTION 


Thermal stresses have become increasingly important in 
engineering practice during recent years. In power genera- 
meomenigher cycle temperatures and use of nuclear fission 
are largely responsible for this trend. This thesis describes 
BEEComputer program for finding temperatures and stresses in 
bodies having axisymmetric geometry and loading. The govern- 
Miemeduations are those of isotropic, uncoupled, quasi-static, 
Bear thermoelasticity. They are discretized by using the 
finite element method. А FORTRAN IV computer program using 
double-precision arithmetic has been written to solve 


problems of the following kinds. 


A. TEMPERATURE PROBLEMS 

The transient temperature vector, evaluated n the nodal 
points, may be obtained for an axisymmetric body with the 
combination of insulated, convection, or constant temperature 
boundary conditions. For the convection thermal boundary 
condition, however, we may have fluid flowing with entry 
temperature prescribed as a linear function of time (RAMP). 
The program can handle up to 15 different ramps, each having 
a different flow velocity, and with discontinuities between 


successive ramps. 


PS TRESS PROBLEMS 
The program will generate load vectors for pressure load- 


ing, centrifugal loading, and axial force. Provision is 
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meee το direct input of one additional load vector. Stresses 


may be found for any combination of these loadings. 


O THERMAL STRESS PROBLEMS 

Thermal stresses may be found for as many as 20 different 
Temperature vectors which may be output of the temperature 
problem or direct input. In short, in this part any combina- 


tion of the temperature and stress problems may be used. 
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Mee | bet beMeENT FORMULATION OF HEAT CONDUCTION 
IN AXISYMMETRIC BCDIES 
A. METHOD OF FORMULATION 
For bodies of revolution under axisymmetric loading the 
mathematical problems presented are two-dimensional. The 


governing equation for non-steady heat conduction is 


КҮТ = ocT, (1) 


<l 


Té Kk is the thermal conductivity, p the density, c the 
specific heat, T the temperature and V the gradient operator. 
The superior dot denotes a time derivative. 


plying Сатеткіп' 5 principle [1] gives 


f N,.V:kVT ау = / oc №. Т ау, (2) 
= ү u 


Mire” the integral is over the volume V of the conducting 
body and N; is a "shape function" used in representing the 
Femperature distribution. If S is the surface of the body 
and n the outward normal to surface, then using Gauss' 


theorem we can write 


τ. 7 E 9T 
E V (N; кугу αν : N; Кт dS (3) 
Since 
τ... τς ay = $ (vN:)*(CkVT)dv 
V E V ; 
+ f N.V- (kVT)dV, (4) 
y 1 


we can combine Eqs. 2, 3 and 4 to get 
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: = eic _ oT 
x pcN; TAN + E (VN: ) E Ni ky 05% (5) 


Each node of the solid region has a separate discretized 
саг equation calculated from Eq. 5 using the appropriate 
Shape function N;. Π СЕО о the volume integrals on the 
πι 5 τ᾿ ΟΤΙ. 5 yields a square coefficient matrix 
in the assembled set of equations. Calculation of these 
Matrices 1s a standard process. Details are given by 
Zienkiewicz [1]. 


The discretized set of equations takes the form 


Toyo. GM 


Ia 


T + 


|< 


Kun there is a term by term correspondence ατα κα.» 
Meal symmetric matrices C and Y represent, respectively, 
the thermal capacitance and thermal admittance. The elements 
mu coyector T are nodal temperatures. The vector v, dis- 
cussed in the following section, depends upon the thermal 
boundary conditions. 

In the present development piecewise constant material 
properties have been assumed, i.e., k,p and c are constant 
within each element, but may vary from element to element. 
Also two-dimensional isoparametric elements are used.  Appli- 


cable equations are summarized in Appendix A. 


m this text a double underline denotes a rectangular 
matrix and single underline denotes a column vector. 
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B. THERMAL BOUNDARY CONDITIONS 
ihermal boundary Conditions affect only those scalar 
equations of Eq. 6 which correspond to boundary nodes. 
Exosdingly, the vector v is sparse. Also, in a single 
Eucblem 1t 1s common to have different thermal boundary con- 
ions On individual portions of the boundary. In what 
Bellows the subvectors of v (distinguished by individual 
superscripts) which correspond to separate boundary condi- 
mons are treated individually. 
1. Insulated 
It is clear that for insulated boundary conditions 
the subvector y UJ αν. τμ hand side of Eq. 6 corres- 
Bondıne to this boundary condition is zero, since = = 0. 
2. Convection 
The heat transfer mechanism occurs in the interface 
Mene solid and fluid. If the fluid temperature is 0 and 
the heat transfer coefficient is h, then equating heat con- 


ducted away from the surface to the efflux from the solid 


fj} gives 


οτι _ 
-k (3), = h(T-0), (7) 


- 


where the subscript S means that the derivative is evaluated 
at the surface. 


Eor constant h: 


ΠΠ τα I 
f N; kon d$ hf NC T) ds. (8) 
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that tollows the fluid temperature 0 is taken to be a 
ecified function of position and time. For purposes of 

EE uwetrzatron; the fluid temperature is specified at a 
discrete number of fluid "nodes." I£ Ө, represents the 

υιό temperature at fluid node j, then the fluid temperature 


Eus the boundary may be represented by 


0 - X N.0,, 9 
18: (9) 


where the N; are one-dimensional forms of the shape functions 
ШЕСЕ for the solid. Substitution in Eq. 8 gives 
pn. k 22 as = 5 y:.(0.-T 10 
: τς = Fag E 13» (10) 


ς 1 ən ; 


Bye the summation extends over the surface nodes and the 


* 
coefficients E are given by 


x 


Assembling the contributions from Eq. 10, the subvector 


ν (2) for the convection boundary condition may be written 


x - Y'e. (12) 


Mme contributions -P Yit; Т ТОП ЕДТ 10 are Included by 


Acne matrix Y (see Eq. 13 below). 
Os tant Temperature 
If 0 represents the constant temperature desired at 
the wetted surface, then we can use the convection boundary 


20 


condition and replace h by a big number (say 10%”). Since h 


ΠΕΙ, large, then for thermal equilibrium the temperature 
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M (he surface will be forced to equal 0. So the subvector 
y O) for the portion corresponding to the constant temperature 
Boundary condition can be obtained from Eq. 12. 

Uo hE aplication Of these boundary conditions in 
INTE le problem, the right-hand side vector v will be com- 
bined from the corresponding subvectors and the finite ele- 


ment discretized equation becomes 


5 


HL VS (13) 


la 


T + 


l< 


+ * 
Кисе у = тҮ + 


[< 


EXACT TIME SOLUTION WITH SPATIAL DISCRETIZATION 
υπ ος Only the solution of Eq. 13 for v = constant 
σι τα at time t = 0. Let de be a particular solution 


(steady state) with ША = 0 so that 


ip Ca | (14) 


For the homogeneous equation 


+ 


T + 


[e 
|< 


пе 0 (15) 


m assumption | = w exp(-At), where w is a vector and À is 
a scalar, yields the form 


+ 


l< 


w = À 


Ia 


W. (16) 


It is apparent that Eq. 16 defines an eigenvalue problem. 
Let A be the spectral matrix and W the modal matrix with 


Normalization according to 


T | (6123) 


I= 
In 
I= 
I 
I 
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mere I is the identity matrix of the same order as C. Now 


let 


EX (AT) b (18) 


Mme Db is a constant vector. Substituting this in Eq. 15 
gives 


+ 


A m LE (19) 


I< 


Шер Ет ЕБ = 


[e 
|= 
|> 


NONE. 19 is satisfied for all b if 


Түу = (19') 


> 


l= 
[< 
l= 
I> 


mhis 1S guaranteed to be satisfied since A and W are 
Spectral and modal matrices for the eigenvalue problem of 
Pope with W normalized according to Eq. 17. 

Foturning to the original problem (Eq. 13), the com- 


mete solution may be written as 
ПЕЕ Е ерла) Ы) (20) 
where 


T = 
=> 


|= 


В, (21) 


and ß is a constant vector. Now B may be found (using 


mae 14) to be 

pU. (22) 
As tituitine this result into Eq. 20 and using the initial 
condition gives the result 

b=Wla- Alu y. (23) 

Meee enerali solution of Eq. 13 may thus be written 


as 
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T 


Т = W(L-exp(-A0) A 3 W v + W exp(-AT)W (24) 


ο Ὁ ος Of the present program non-zero compo- 
ES Of vector v are to be specified as piecewise linear 
me ons of time, During each segment of time history of v 
Euunalytrical solution of Eq. 13 is possible in the form of 
a particular solution plus a complementary solution such as 
E0. At each node the corresponding time variation of 
temperature wlll consist of a linear part contributed by the 
particular STET and a sum of n terms representing the 
uudementary part. Each of these n terms decays exponen- 
tially with a separate time constant. In principle it is a 
Isi Cíorward process to find each particular solution and 
accompanying complementary solution. 

Contemplated problems may typically have from 10 to 
40 segments required to represent the piecewise linear 
variation of v. The number of body nodes n will be of the 
order of 100 or more. In view of the number of particular 
solutions required, each accompanied by an individual comple- 
ШОО ату solution of the form given by Eq. 20, it was concluded 
that a numerical solution of Eq. 13 would be considerably 
Were ceconomical than an analytic one such as that given by 


τα. 24. 


D. TIME INTEGRATION 

In this section the relative merits of the Runge-Kutta 
and trapezoidal methods of time integration are discussed. 
Since either of these methods will give an exact result if 


tle solution is a linear function of time, investigation 
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. υπίιπεα to performance on a single scalar equation 


ISO (25) 


whose solution y = n exp(-At) is of the same form as the 
components of the complementary solution (Eq. 20 ). 
M kunge- kutta Method 
A method introduced by Runge and subsequently elabo- 
wed by Heun and Kutta [3] is widely used for the numerical 
ton of first order ordinary differential equations. 
MS algorithm prescribes a Sequence of calculations for 


determining the ordinate У; at Cime T: 


= ¢.+At in terms 
+] 1+1 ti τ 


of m and values of y at intermediate and end points of the 
miterval At. The fourth-order form, which requires four 


evaluations of y, gives for Eq. 25 the result 





Yi+] (A49)? _ Our)” „ (OAT) 
КИСИ τπτ со 


ehrt hand side of Eq. 26 represents the first five terms 
of the Taylor expansion of the exact solution 

Y = exp(-AAt)) so we may conclude that the relative 
Sagem in cach time step is less than modulus of the next 
term: (ААт)?/51. 

In addition to providing the apparent prospect for 
high precision indicated by this error bound, the Runge-Kutta 
method also permits changes of time increment during the 
integration process without requiring additional computa- 
tionally expensive matrix decompositions. The attractiveness 
of these two features μοι a thorough exploration of the 


potential of this method for the present application. The 
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meecontined to performance on a single scalar equation 
y + Ay = 0 (25) 


whose solution y = Y exp(-At) 1s of the same form as the 
components of the complementary solution (Eq. 20 ). 
is kunge-Kutta Method 
A method introduced by Runge and subsequently elabo- 
mated by Heun and Kutta [3] is widely used for the numerical 
Solution of first order ordinary differential equations. 
This algorithm prescribes a sequence of calculations for 


at time T. 


Setermining the ordinate DE е 


= . + 1 
-ῃ DE AT Jj) terms 
of y, and values of y at intermediate and end points of the 
u val Ντ. The fourth-order form, which requires four 


evaluations of y, gives for Eq. 25 the result 





DENS а) EET 
EN ou 39s C9 


e Ticht hand Side of Eq. 26 represents the first five terms 
of the Taylor expansion of the exact solution 
Yiri = exp(-AAt)) so we may conclude that the relative 
error in each time step is less than modulus of the next 
term: PULSES 

"ΠΤΙ. τυ providing the apparent prospect for 
SN precision indicated by this error bound, the Runge-Kutta 
method also permits changes of time increment during the 
1ntegration process without requiring additional computa- 
tionally expensive matrix decompositions. The attractiveness 
of these two features ee a thorough exploration of the 


potential of this method for the present application. The 
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disqualifying defect which emerged after studying a number 
of examples is readily appreciated from examination of 
Table I. For values of AAt less than 0.5 it is apparent 
Mat Runge-Kutta scheme affords acceptable engineering 
accuracy. However, when the method is applied to solution 
KE 15 we must deal with a number of A's equal to the 
liber of nodes (see Eq. 20 ). This number may be greater 
than 100 and the ratio of the largest X to the smallest may 
Eu cexceed 1000. Although the solution is dominated by 
the contributions of the eigenvectors corresponding to the 
smaller As, it is clear that the solution will be unstable 
if the largest ААт exceeds about 2.7. Because an unaccept- 
Hub small Av 1s required in typical problems, the Runge- 
Kutta method was rejected. 

2. Trapezoidal Method 


The trapezoidal method estimates yi,] from the formula 


_ DC | 
РУ о (у τ)’ (27) 


Substituting for y; and m ivemsta. 25 and rearranging 


gives 


m: E 
laz AAT — (28) 
yi ο. AAT 


ΙΙ expansion of the right-hand side provides an error 
bound (per step): 
5 
(AAT) /12. 
From the point of view of the size of AAT this method has no 


stability limit, but has slow attenuation with alteration 
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TABLE t 
ATTENUATION FACTOR CCMPARISON 


Fourth Order Runge-Kutta Algorithm 


AAT η; exp (-AAT) 
(Runge-Kutta) (Exact) 
.0001 . 9999 . 9999 
‚001 ‚9990 ‚9990 
.01 . 9901 . 9901 
51 .9048 .9048 
32 S157 291557 
215 .6068 .6065 
1:0 29750 . 3679 
20 105555 555 
2.5 .6484 0821 
3.0 123750 .0498 
4.0 5,0000 204183 
5 ῳ 1102535535 .0003 
ο 0 291.0000 .0000 
20.0 E519 555 .0000 
50.0 240784.3333 .0000 
100.0 4004901.0000 .0000 
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Econ for large AAT. Table II shows this behavior. 

Since a wide usable range of AAT is essential and 
the stability of trapezoidal integration is guaranteed, 
Bhus method is S cen for the present program. 


Applying the trapezoidal algorithm to Eq. 15 yields 


A girl = G Ti + AT (yitl + vi) (29) 
= «απο --- -- 2 == 3 
where 
A = С + AT Wiss 
— — 2 ws 
S TET 
G C > Na: 


and the Ber үрт» denote evaluation at discrete time 

M ruas At. If m is the order of the capacitance and 
τις: matrices, C and Y, then once a certain step size 
Ат 15 chosen, it requires m?/3 Operations to perform the 
Meeded triangular decomposition of A. Thus, ο m, a 
change of step size At becomes costly from the point of view 
of computer time. Accordingly, in the present program only 
ишетте step Size is used throughout each problem. 

Also, for assuring sufficiently rapid attenuation of the 
components corresponding to the large AAT, the following 
correction is utilized. 

Sons Correction 

Irons proposed a scheme [4] for augmenting the atten- 
uation of the contributions of those eigenvectors for which 


AAT 1S large. Define 
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TABLE TI 
ATTENUATION FACTOR COMPARISON 


Irapezoidal Integration 


us Ta ELLE 

(Trapezoidal) (Exact) 
.0001 . 9999 . 9999 
001 . 9990 . 9990 
ο] . 9901 ‚9901 
‚1 ‚9048 ‚9048 
22 8182 „81867 
6 2259] .7408 
5 . 6000 .6065 
0 15935 . 3679 
0 . 0000 1.558 
5 ΠΠ ‚0821 
τι 22000 .0498 
.0 2035555 20185 
.0 -.6000 .0003 
. Ü -.6667 .0000 
¿0 = 8182 .0000 
‚0 22925] .0000 
. Ü -.9608 .0000 
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* 
y ade 25 р SL Ὃν, 


1 iti? (30) 


where y and y- 


j+] are obtained from Y4 -1 by trapezoidal 


Mitegration. 
tie program presented in Appendix B, Eq. 30 is used 
Mater every 10 steps of time integration. Table III shows 


the resulting modifications. 


E ESTIMATION OF EXTREME EIGENVALUES 
Euadtrc results for one-dimensional heat conduction 


ive, for an eigenvalue 
3 3 


л > | (51) 


where d is the distance between points of extreme tempera- 
Ше апа zero temperature. 
If we use this to estimate the smallest A in Cylindrical 
coordinates, two modifications are recommended. 
resume that the point of zero temperature is in 
the fluid at a distance from the wall equal to k/h, 
where h is the surface heat transfer coefficient. 
2. If there are two approximately orthogonal paths for 
heat flow from the (single) maximum temperature point, 


then replace PT in the above formula by 


2, "1. ...; Y | (32) 
d a2. d 
min. max. 





PFOormecotimatine the largest À, the surface heat transfer 


ο στ πο μας πο Significant effect. We may continue to 
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TABLE III 


ο δν ΤΟΝΙ CORRECTION 
TIER IO STEPS OF INTEGRATION 


AAT рало) 710/70 Y" 10/Yo 
Exact Trapezoidal Corrected 

0.00010 ООО 0799900 0.99900 
0.00020 799300 0.99800 0U ου 
0.00100 = 21005 0.99005 0.99005 
0.01000 .90484 0.90484 0.90486 
0.10000 290.798 0 56757 0.36849 
0.20000 .15534 0.13443 0. 15579 
0.30000 ο 0.04860 0.04978 
0.50000 ‚00674 0.00605 0.00645 
1.00000 .00005 0.00002 0.00002 
2.00000 .00000 0.00000 0.00000 
4.00000 .00000 0.00002 -0.00001 
8.00000 .00000 0.00605 -0.00040 
10.00000 . 00000 0.01734 -0.00072 
20.00000 .00000 0.15445 -0.00136 
50.00000 .00000 0.44914 0.00077 
100.00000 .00000 0267028 0200027 
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use the same formula for =; but now consider only the 
d 
smallest element and take 


d Pzyensth orzsmallest side 


min A 


(32') 
d ici enon hanes es side 


max 4 


ENonpargsson of estimates based on Eq. 31, 32, 32' with 
E Exact Solution Of the eigenvalue problem has been carried 
out for several examples. Based on these comparisons, it is 
believed that these estimates are sufficiently accurate for 
Emoosing a time step and estimating the time of occurrence 


Oi the extreme stresses. 


28 





ITI. ONE-DIMENSIONAL HEAT CONDUCTION 


For comparison of numerical (time) integration methods, 
Studies of one-dimensional heat conduction were made. In 
ЕЕ section numerical results for trapezoidal time integra- 


ШО и<їпр Irons' correction are compared with the exact 


transient temperature solution. 

mm ilet ar flat slab of thickness L with conductivity k, 
EN Pp, specific heat c, zero initial temperature, one 
mee insulated and the other in contact with fluid at tem- 
perature T, ЕОС еза асе пеат transfer is h. The 


ЧЕ transient heat conduction solution is available [5] as 






o sinuiL Cosu;x “us x 
T=T 1-2 Σ —— E (55) 
f ve | u;L*SinuiL CosuiL 
» 
mS та 


— "Temperature Te 


= h_ -_ - — — 


Figure 1. Slab with One Face Insulated and Another 
ime ontact vith Fluid., 


23 





Fire k/pc is the diffusivity and the u; are the solutions 


of the transcendental equation 


zonL A 
Tan uL = k NIE 


бз 
lene Limite element comparisons we subdivide the 
Eur tance L into m one-dimensional 3-noded elements and use 
the corresponding isoparametric shape functions. (The work- 
ing equations, shape functions and the element capacitance 
and admittance matrices are given in Appendix A.) 
Mathe course of this investigation separate computer 
programs were written to evaluate nodal transient tempera- 
tures using the following methods: 


(а) exact transient temperature solution, Eq. 33; 


(Б) “exact time solution with spatial discretization 
section 11,07 Eq..24); 


ο ре Kutta time integration (section II.D, 
part a); i 


(d) trapezoidal time integration (section II.D, 
part b); 


ΙΙ. |ozommdaletrmeezntesration with Irons' 
connection. (section II.D, part c). 


Por an initial step change of fluid temperature from 
Zero to 1, transient temperatures have been found. For these 
comparisons the parameters (in consistent units) were taken 


to be: 
L = 8, р = 25, k = 8, c = 5 and h = 5 


The distance L was subdivided into m = 3 elements. For the 
present purpose, comparison is confined to the exact solu- 


tion (item (a)), and the finally adopted system (item (e)). 
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In Table IV, temperatures at the two faces (x = 0, x = 8) 
meee compared for various times. The trapezoidal integration 
Mas Deen performed using the constant time increment unity 
рша спе Irons' correction is applied after every 10 incre- 
ments. 

It is believed that Table IV demonstrates that the 
merical integration method gives adequate accuracy for 


engineering applications. 


TABLE IV 


COMPARISON OF ONE-DIMENSIONAL TRANSIENT TEMPERATURES 


Exact 


Trapezoidal 
with Irons' 


Exact 


Trapezoidal 
with Irons' 





ons" correction is used after every 10 steps of trapezoidal 
integration. 


ο 





οσο PROBLEM 


for bodiles of revolution deformed symmetrically under 
axisymmetric loading, the stress distribution is two- 
mensional. Since deformation is symmetric about the axis 
EN cVUolutron, cylindrical coordinates (R,Z2,0) are used. 
Eum llows that the stress components are independent of the 
mere O and all derivatives with respect to 6 are zero. 
"io the components of shearing stress tno and 779 vanish 
puaccount of the symmetry. But since any radial displace- 
ment induces a strain ερ In tie circumferential direction, 
mE non zero component of strain and the three in-plane 
components C Ep > YRzJ> eomplete the state of strain at a 
femme in any axisymmetric situation. Hence the state of 
stress for an axisymmetric body under axisymmetric loading 


is given by 


DE MU EE: (55) 


Z 
its chapter the stiffness matrix of an axisymmetric 
body and the thermal, pressure, and centrifugal load vectors 
ero mutated and, finally, evaluation of stresses at a 
eomte is discussed. The treatment closely follows that of 
Zienkiewicz [1] and this reference should be consulted for 


Further details. 
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КОО STIFFNESS MATRIX 

Ihe elements used are bodies of revolution (about the 
Mois). Ест analysis it is sufficient to describe the 
wes section in the R,Z plane. In Fig. 2 such an element 


the local £,n coordinates are shown. 





Figure 2. Quadrilateral Element Representation. 


If u and w are the displacement components at a point in the 
directions of R and Z respectively, then these displacement 
components may be defined in terms of the nodal displacements 


by the appropriate isoparametric shape functions as 
8 8 
q 5 миг, w= X N.w. , (36) 


where N; , a function of Ë and n , is the shape function for 
Clement node i and Ui,Wi, are Масите са атериасешеп е compo- 
nents. The strain-displacement relations [6) can now be used 


to obtain the components of the strain vector. Thus 
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25 


[5ο 


δ; (37) 


Mere Bis the standard rectangular strain-displacement 
Meer ix of any finite element formulation, a function of 
ES local coordinates š and m, and ë is the vector of nodal 
placement. (See Appendix A, part 3, where the applicable 
formulas and useful equations are summarized). 

EN Ho olasstucyty matrix for an isotropic material is 


then the stress vector c at a point is given by 


oa 


| 


ε. (28) 


Now, by evaluation of the total strain energy in the element, 


WS element stiffness matrix can readily be obtained as 


ESTEBE DBEa, (39) 


Meere the integration extends over the volume of the 
ement. 

In the present program the upper triangle of each ele- 
Ment stiffness matrix is evaluated by numerical integration 
using four Gauss points within the range of & and of n [1]. 
Me lement contribution is immediately placed in the system 


San ness matrix, which is stored in banded form. 


Eee THERMAL LOAD VECTOR 
If we denote T as the difference between local tempera- 
ture and reference temperature, then the thermal strain £5 


1S given as 


Da, (40) 
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where 


πο 


and a is the coefficient of thermal expansion. 


The thermal load vector jg 15 given by Zienkiewicz [1] 
as 
M ES Dee dV (41) 
—T = == ° 


ШОШ che point of view of the numerical evaluation it is 


Mieeresting to note, however, that the product D oe in 


Eq. 41 will reduce to 


O 
tri 
E 


D to  ]1-2v ы» > 





where E is the modulus of elasticity and v is Poisson's 


ratio. Thus 





Siu ГО T 
Fr wen I N (43) 
n 
where T = X N; T; and n 1s the number of nodes in each 
1Ξ1 


element. 

In the attached computer program in Appendix B the advan- 
mee Oot the simplicity of Eq. 42 has been utilized. Also, 
КО е Бас some zero components, in the process of multi- 
plication of в! Αν Е асо топ ов еле appropriate 
Шоле гето components of each column of the Pi hasi been iper- 
formed. Finally, Eq. 43 has been integrated with four Gauss 


points. 
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O PRESSURE LOAD VECTOR 

Consider a quadrilateral element as in Fig. 3 on the 
boundary of the axisymmetric cross-section where constant 
mmessure P 1s applied. The infinitesimal force dF due to 
mae normal pressure acting on the inner infinitesimal cir- 


@amrerential surface dS is 
dF = PdS = 2TRP dí , (44) 


where d£ is the infinitesimal length along the side of 


quadrilateral. 





Z 


Етрпте 5. Boundary Element Under Pressure. 


кет ЕВ апа F; be the components of the pressure force in the 


Ed Z; directions respectively, then 


dF, = dF Cos 6, dE, dl SIN CL, (45) 
where Sin 0 = = and Cos 9 = 92 . Therefore 
ЧЕР τρια, dE, - -(2m RP)dR. (46) 
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Since the total work W done by the normal force F is equal 


ο the Sum of the work done in R and Z directions, 


W fu dF, + fw dF 


2 


РОЧ Ев мав). | (47) 


EE by using the appropriate shape functions, each compo- 
mm Of Eq. 47 may be defined in terms of the nodal values. 


Since 


T 


М = (6°) Е5, (48) 


where nc NNNM c lem O nthpressure load vector 


ESI contributed by node j is explicitly given as 


e ON; 
Er 1 Ll 
ЕС = 21р f (ЕМ. А.) N; dE, (49) 
р |е, -1 aN. 
Jp ER R. 


where N; in Eq. 49 is the appropriate shape function for 


node j. 


Meee CENTRIFUGAL LOAD VECTOR 

БО ст Каргаш Fig, 2eand assume that the body is rotat- 
ing about the Z axis with eu la Speed еп the centri- 
Lugal force per unit volume at a point distant R from the Z 
axis will be on? R, where p is the density of the material. 


The work done in this case is 


Кс л. (50) 


V 


o 





Bor constant ES the corresponding element centritugal load 
vector is readily obtained by evaluation of the components 
of the integral in Eq. 50 in terms of the nodal variables, 
me. , 

j 


B - 2702" J 


2 
j. (IN;R.)^ det J N; d£dn, (51) 


Z 
Iu 
Sere det J is the ΠΠ of the Jacobian coordinate 


transformation matrix (see Appendix A). 


E. STRUCTURAL BOUNDARY CONDITIONS 
The structural boundary conditions implemented in the 
program are: 
(a) one or more nodes prevented from moving axially; 
(b) one or more nodes Prevented from moving асаа ту 


(c) the right-hand end cross-sectional plane remains 


The computer program presented in Appendix B has the 
capability of handling any combination of the above men- 
tioned structural boundary conditions. For boundary condi- 
Pons of the types (a) and (b), Simply multiplying the 
corresponding diagonal Component of the stiffness matrix by 
1020 gives zero displacement [ο]. TOT practical purposes). 

For the boundary condition of type (c) both ends are initially 
fixed axially for all solutions. An additional ΟΠ S 
obtained for unit axial displacement of one end. еа тат 
force is evaluated for each solution. Superposition is pers 


formed by adding the displacement vectors for the given 
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ings (thermal plus mechanical), plus an appropriate 
Euetron of the vector found for unit axial displacement. 
NUS fraction is chosen so that the resultant axial load has 


Mm specified value. 


E SYSTEM EQUATION SOLVER 

mee the desired structural boundary conditions are 
ed, then the problem is to find the nodal displacement 
EN or, corresponding to a given number of load vectors. 


We have 


IR 


er (52) 


Encre K is the system E ne e matrix in banded form and F 
1s a load vector. In the present computer program a single 
Mes ky decomposition is performed on K. Then, by a process 
i orward and back substitution, each load vector is 


ὅσα by the corresponding displacement vector. 


PPI PINCIPLE OF SUPERPOSITION 

Upon the evaluation of the displacement vectors due to 
EME various types of loading, the principle of superposition 
can be applied on the displacement vectors. On each thermal 
displacement vector the displacement due to any other type 
of loading is superimposed and, as the result, the number of 
displacement vectors is reduced to the number of thermal load 


vectors. 


E SIRESS EVALUATION 
EXonmcEtemedgrsplacement vector ó, the displacement 


СЕВО ol each element δ΄ may be obtained easily. Then the 
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: ° e . 
wene ponding element strain vector e` at any point can be 


found from 


e 


zt oc mod (53) 


[το 


Bay, the corresponding element stress vector σ΄ 15 


obtained by 


of " < Ic (54) 


Є 
πο 


Io 


Since normally the stresses on the inner and outer surfaces 
EN rsymmetrrizc bodies are desired, in the computer program 
presented in the Appendix B provision has been made to cal- 
culate the stresses at the two Gauss points corresponding 


to Ë = + on The inner and outer boundaries of each 





ν 2 
element (where n = #1). 


Upon the evaluation of the stresses at each point the 
mean stress and the octahedral shearing stress [7] are 
lated. The program gives as output the extreme values 
"δες stresses, the R and Z coordinates of the corres- 


Ending points, and the times of occurrence. 
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ΘΙ AD TRANSIENT STRESSES 


In this section a one-dimensional comparison of stresses 
1s made between exact and finite element results. The tran- 
Eu temperature problem 15 the one previously described in 
Section III. 

hé slab edges are free to translate in the plane of 
Elle Slab, but are nC vented fTOM rotating, the exact solu- 
tion for thermal stress [9] is 


M ENS В 
Е те ως πο (55) 


where T is the local temperature (Eq. 33), and 


Tave ute avéerape temperature in the slab 


If we choose E = 2, a = .50, and v = 0 (all in consistent 
nau then the maximum stress obtained from the exact solu- 
tlon 1S 


ee = -,477786 


mau It occurs at х = 8, т = 73. 
Using the finite element technique with trapezoidal time 
integration and Irons' correction every 10 steps, the maximum 


E ress is found to be 


ο ος. = -.477778 


and it also occurs at x = 8, t = 73,.as before. 
P—iomoncenvedsthiateche method chosen gives excellent 


results. 
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Vie leon PROBLEMS 


Program integrity and accuracy have been verified by 
Serving a number of test problems. Since stresses, whose 
evaluation depends upon derivatives of displacements, are 
known to be less accurate than temperatures, comparisons 
ES Exact results are confined to stresses. Individual 


Broblens are described below. 


ESSSIHICK CYLINDER 
monomer a thick cylinder with inside radius 30 inches 


and Outside radius 50 inches and the following material 


[moeperties. 
Modulus of elasticity E = 28.9 x Tob Psi 
Poisson's ratio v = .28 . = 


7.22 х 100 1/9 


Meer ticient of thermal expansion & 


AE 8 Btu 
Thermal conductivity k = 28. a 
Density p = 489. Lbm/ ft 
Specific heat TT SEU 
p à ШОШ К 


An arbitrary length of 25 ones has been selected for the 
sender and it has been subdivided into two different 

ου πει representations as in Figs. 4 and 5. The plane-end 
boundary condition with zero axial force is used. The 
Stresses tor various types of loading are compared with the 


Epmesponding exact analytic solutions as described below. 
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co Radial Elements Representation 
ol Tick Cylinder. 


e Radial Element Representation 
Thick Cylinder. 


1. Thermal Loading 


Бек елт Куат таш оп о Тїе temperature T = ZOR 
u (C rosses obtained by the finite element method for the 
above representations are compared with the exact analytic 


Solution in Fig. 6. The 1 oblea clearly 15 


9 


RZ 


zero and the one obtained by the program was 10. The 


racy Of the Other results is clearly satisfactory. 
ES Sure Loading 
AO pressure Of 1000 psi. acts on the inner 


sur ace of the cylinder. Again, T is zero and the program 


БОК 


RZ 


gives 10 In Fig. 7 the other stresses induced by this 
uniform pressure are compared with the exact solutions. 
Here also the accuracy of the results, even with two radial 


elements, is adequate. 
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ЛСТ 
A 2 ELEMENTS. 


NO 


og > 





e 


FIG. 6 
THICK CYLINDER UNDER THERMAL 
LOADING, T= 20R 


44 





PS! 


2000 So 


1600 


1200 


—— EXACT. 
800 A 2 ELEMENTS 
O 5 ELEMENTS 





400 


-400 


- 800 





Pu ccc G 35 E 


Fig.7 THICK CYLINDER UNDER INTERNAL PRESSURE, 


45 





E Centrifugal Loading 


The speed of rotation has been assumed to be 500 
Evolutions per minute. The stresses obtained by the pro- 
mm are compared with the analytic solution in Fig. 8. In 
iS case also the results obtained by the program, even 
meen Only two radial elements, are very close to the exact 


Solution. 


fee HOLLOW SPHERE 

Mec onsider a hollow sphere with the same material 
Meoperties as in the thick cylinder with inside spherical 
radius 30 inches and outside 50 inches. The loading is 
thermal with T = 20 ‘(spherical radius). Symmetry permits 
mening Only half of the sphere for the computer analysis. 
Er elements representation is given in Fig. 9. Since the 
program gives stresses in cylindrical coordinates, these 
have been transformed to the spherical coordinates for 
Comparison with the exact solution in Fig. 10. The accuracy 


Of the results is noteworthy. 


pees JHERMAL STRESSES IN NOZZLE 

This problem concerns thermal stresses near the inter- 
deton Of a cylindrical pipe and the spherical vessel. 
Fig. 11 gives the cross-section of the structure which is 
be analysed. The material properties are: 


ANM S Eb» 


tt 


Modulus of elasticity E 


Poisson's ratio у = ,30 - 


6 


7 б х 10. 1/Ε5 


Thermal expansion coefficient a 
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KSI. 


4.5 


3.5 


2.5 


L5 


=5 





— EXACI 
eee se LEME 5. 
ος ELEMENTS: 


FIG, 8 


ROTATING CYLINDER 


— F D MH — —ÓÓ——Ó— ο 
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Eigure 9, 


Hollow Semi-sphere 32 Elements Representation. 


48 





60 


7O 


60 


50 


40 


30 


20 


IO 


ΙΟ 


EO 


30 


40 


50 


Kol 





307732 34 


Figure 10. 


36 38 340 42 44 46 48 SO SPHERE 
RADIUS. 
IN. 


Thermal Stresses in Hollow Sphere 
I 2 d spHerrcaleradius). 


49 





Z, ONYT, 


S3Ow3ell$ N33^.138 
NOILv2ay Vm 3O 


„2, 3004805 


81081 
Ol ΓΈ 





3790109 Y 8208! 


οι ονν 37049 ὃ 5' 01 
LNZONVL ЭМ LHOIVALS 





AGL AINOWD FIZZON 
ДО 





so 





Thermal conductivity ΙΤ SDturftt.hr.?F 
Density p = 530.5 Lbn/ £t? 
Specific heat CoL Dux D bm? 


The loading results from the thermal transient in the fluid 
Koli Ined in the nozzle. This fluid is in contact with the 
Eure on surface "I" (Fig. 11) and has the entry temper- 
îre time variation as in Fig. 1Z. The fluid in the sphere, 
Mis in Contact with the structure on surface "2" (Fig. 
Mas the constant temperature 478°F. At t = 0 structure 
EN uniform temperature of 478°F and is stress free. 
Mereni1or Surface of the structure is ον, The flow 
velocity past surface 1 is 8.5 ft/sec (inward) and 
Biere is no flow past surface 2. Τις οσο ο ек 
Bramster coefficients are 1395 and 2910 τ ГОТ 
ШО Гасе< "1" and "2" respectively. 

Bor che structural boundary condition it is assumed 
mat the nodes on the left end cross-section of Fig. 13 are 
Beevented from any axial displacement. 


The maximum thermal stresses obtained by the program 


pecur in element 14 of Fig. 13 as follows: 


Maximum mean stress = 18.81 ksi. at time = 4 seconds; 
Maximum octahedral shearing stress = 11.5 ksi. at time 
= 18 seconds. 


These results appear to be reasonable, but no suitable com- 


parison solution is available. 


Sl 
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(I CONCEUSTONS AND RECOMMENDATIONS 


fae CONCLUSIONS 

Pecomputer program has been developed for the solution 
Of axisymmetric transient heat conduction and thermal stress 
problems. The system will accommodate a wide variety of 
Beometric arrangements, thermal and structural boundary con- 
ditions, and mechanical loadings. 

Although double-precision arithmetic is employed through- 
out, efficient algorithms for the manipulation and storage 
of large symmetric banded matrices allow in-core solutions 
with modest time requirements. 

πο ουσ τατις isoparametric elements used allow accurate 
representation of curvilinear boundaries and the stress 
field. Examples presented show that a small number of ele- 
ments is generally sufficient to determine stresses with 


pood precision. 


B. RECOMMENDATIONS 

Incorporation of several additional features would 
ο Παπ icantly increase the program capabilities. The follow- 
ing extensions are recommended. 

Есета thermal and elastic properties have been 
assumed constant within each element. Since such properties 
are generally temperature dependent, provisions should be 
made for uU E daring them during both temperature 


ana. stress solutions. 
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2. The surface heat transfer coefficient, taken as 
u ant in the program, is a function of temperature and 
flow velocity. Provision should be made to include these 
effects. 

Pee ine program presently starts every temperature 
10η with constant initial solid temperature. Provision 
for an externally specified initial temperature vector should 
Бе included. 

ШИ еа ое Чап ту бї temperature and stress results 
mem@erated Dy the program is currently presented as digital 
EUtout. Graphical output in the form of two-dimensional 
contour plots of temperatures and stresses should be 


mmovided. 
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APPENDIX A 
APPLICABLE FORMULAS AND EQUATIONS 


PART 1 


(a): Two dimensional 8-noded (parabolic) 
M oparametpreoshnape functions. 


Corner nodes: 


N; = + "mon ΝΗ; 


Mid nodes: 


R Ee: - | 
5i: 7 0, N. = 2 (1 ΠΛ 
= m 2 
n; = 0, Noo oe FE) η΄) 
where 
EES QM ne. 
E Temperature at a point in terms of the 
nodal temperatures. 
8 
T= EX М.Т. 
. 1 1 
1=1 
C Coordinates at a point in terms of the 
Modal coordinates. 
8 
R= E N.R 
ΠΠ 
1-1 
8 
Z= È N. Z; 
1=1 
(d): The Jacobian coordinate transformation 
matrix. 
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(e): 


(Г): 


(g): 


PART 2 


fa): 


(b): 


97 ƏR 
38 dE 
BE 
= 92 IR 
on on 


Element of the finite element capacitance 
Matrix., 


e --- 
E = pc S N; N: dV 


Element of the finite element admittance matrix. 


e = V Oy 
E k J VN; VN, dV 


Elemental volume. 


ue qetu ας an 


One-dimensional 3-noded (parabolic) iso- 
Par metric shape functions: 


End nodes N; = 5 ES а= Eo) 
_ 2 

Mid node N; = (1 FES) 

where, again, BEE 


One-dimensional capacitance and admittance 
matrices. 


4 2-1 70728 1 DEOS 
ns k 

E M RM КЕ 0 

-1 2 4 l -8 7 0 0 


where & is the length of the element. 
The element Vevector for the one-dimensional case. 


€ T 


Uu Ξ <hT £ 0 0» 
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PART 3 


(a): 


(b): 


(С): 


(4): 


Strain- displacement relations. 


Fu 21 


ee s pi s i 


ON. ON. 


с асч cn _ 1 1 
Dr ma nn 





where N; are the same isoparametric shape 
tunc tion as in Part I, (a). 


"ΓΡ ΠΕΡΙ. 


ON, ' ΟΥ; ON, 
με... ч 
IN 
Z 
0 — ОСИУ. 0 
Е i 
= N 
2 
0 B NEM. us 0 
9NJ oN. IN, ΟΝ 
ἫΝ V n: 





The element displacement vector. 


1 l-v 1-۷ 0 
Б E(1-v) vO V 
D (1+v) (1-2v) Ι1-ν 1 l-v : 
i V V 
l-v l-» 1 0 





COMPUTER PROGRAM 
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RETURN 
END 





ΤΠ τ 


USER'S MANUAL 


INCAS. section the procedure for the use of the present 
water program is described. It is intended that a person 
ET minimum familiarity with the details of the finite 
Ei uento method and computer programming be able to use this 
program. 

Ihe steps below are in order and the user is advised to 
aou them carefully. 

For finding the thermal stresses in an axisymmetric body 
.. symmetric loading, the user has the choice of using 
Ne: ο τις or British system of units. The program 
will handle both systems and the necessary conversions are 
made automatically within the program. ООО ТЕСШ тт 
num each system must be consistent and they should be as 
listed below in Table V, 

lector opreparation of the data input for a given axi- 
Symmetric geometry with prescribed thermal and structural 
neusdary conditions we go through the details with a simple 


Example . 


seep 1: 
Mente longitudinal cross-section of the body to scale. 
ШОС пе cylindrical coordinates R and Z, with the 


©rieim Of the Z axis on the most left-hand point of the 
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GABLE. ҮМ 
UNITS FOR INPUT DATA 


Does can input card specifies whether British or Metric 
@ata 1S being used. 





m — 


Quantity rss aen qt. oy Sen Metric Syst. 








Beordinates ТИ cm 

Time sec. c c 
Temperature °р EC 

Velocity Iu Sec. m/sec. 

Axial once! 1bf kg 

Pressure D kg/cm? 

Density | 1bn/ft? gn/ cm? 

ООС Гтс heat Btu/lbm°F Cal KOL 

OM Of elasticity е kg/mm“ 

Bos of thermal exp. ye 17C 

Thermal conductivity D tuuc ευ οκ c j Sec mem oG 
τ transfer coef. BO nor cCal/sec. cm. °C 
Load vector IDE kg 

Rotational speed Revolutions/min. Revolutions/min. 


“ibe ш μπα force and ibm 1S pound mass, 





body. The entire cross section must lie in the first quad- 


maint Of the coordinate plane. 


R 
Io 
© 
О. 
E Эі ΙΟ 15 2 
wre I14. Longitudinal Cross Section. 


ΠΤΙ, ὁ: 

Subdivide the cross section into a maximum of 40 quadri- 
MEE elements. This subdivision is extremely important 
and we should provide smaller elements for the parts of the 
ross section Where we expect the largest stresses or tem- 
ιτ gradient. If the body is made from several mater- 
lals, elements should be chosen so that each element contains 
only one material. Any two adjacent elements must share one 
ο "τσι side of the quadrilateral. Number the elements, 


ИШЕ Шр from 1, in any arbitrary manner (see Fig. 15). 


op ὁ: 


e oli noded elements are used in the program, iden- 
es ese modal points around the boundary of each element. 


ο N les will be at the corners and the other four will be 
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at the mid-points of the sides. Number sequentially these 
nodes starting from 1 and increasing in the direction where 
the number of elements is the least. (The numbers thus 
πυρά are called global node numbers.) For ekayıty of 
this step assume a cross-section as in Figure 15 such that 
the maximum number of elements in one direction is less than 
the maximum number of clements in other direction. (In 

Fig. 15 we have maxima of 2 and 3 elements in R and Z τος. 
tions respectively.) Thus we number the nodes beginning in 
the R direction. See Fig. 16. This method will give the 
minimum band width of the stiffness matrix and will save 


computer execution time. 





Figure 15. Subdivision Inot Elements. 


seen 4, 
At this point we are ready to prepare the first data 


cand Input quantities arc: 
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ΠΕ 0. “Element and Node Numbering. 


NET I total numbertfolí elements 

NNT Total number of nodes 

NCN Total number of corner nodes 

N@FN Total number of mid-Side nodes not lying on the 


Erna lacht line Joining thel corner nodes (number 
Sl E Off’ nodes) 


NMAT Number of different materials 


NPROB Number of problems to be solved 


For the present we take NPRQB - 1. 


Prepare a Single card that reads 
NET, NNT, NCN, N@FN, NMAT, NPR@0B 


Berne format (815). For our example, if all the elements 
BD от the same material this data card reads: 
ee pt 


d Or er MIRAD “ον on ` =. Aa edi rese Ó ` Са eg’ u «> 


00000000°00000 0500000000000000055600000006C0(62006000066GC0006551 
123456891911 213i£15 i6 12 16 189 20 21 2223 24125 26 27 28 23 30 31.32 32 94 35 35 32 33 39 40 41 42 45 4e 45 65 47 48 9 50.51 92 5354599857 5845 δ) 0 
9 9 4 y 4 4 ¢4 + τα ο тор чес ч Бага ч 9 чт ο «1 1.11 4 


4 1 08 4« 9 4 @ # 4 ο 4« Φ. «49: 4 5 
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КОШО е Card, or any succeeding one, if an input quantity 
(e.g., N@FN) is zero, the corresponding field may have a O 


or be left blank, unless otherwise stated. 


ошер ο: 

For each corner node a card should be punched which 
Eds: The node number I, the R coordinate and the Z coor- 
fame Of the node. The format is (15, 2F10.5). 

it 1S not necessary to sort these cards in the order of 
increasing or decreasing corner node numbers, they can be 
Mito cether in any arbitrary order. A typical card is 
illustrated in step 7. There must be as many punched cards 


as the number of corner nodes (NCN). 


Step 6: 

In this step we read in the connectivity array, LE, 
for each element we prepare one card which gives the global 
node numbers and the material identification number. 

For each element, start from the lower left-hand node 
and move in the counterclockwise direction within that 
element and punch the global node numbers in order. The 
tamat 15 (915) where the last I5 is the material identifi- 
cation number. 

It is very important that the cards prepared for this 
ου ο ους together in order of elements, 1.е., тһе first 
Sl oT element 1, the second card for element 2, etc. 

For ease of sorting the connectivity cards the element 
number may be punched after column 50. As an example, for 


element 2 of Fig. 16 the connectivity card should read: 
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ro 


0000000000000000000000000000600002000500 poen pon GOD { 
Г : 3 : : : Í : q n ae τος ο IT 
D QNM EN SENDE ΝΥ S S U p iq t $ 


Memes. identifies the type of material in this element and 


Киа οἱ υπ 55 stands for the element number. 


Бер 7: 

ШО сас side of every element is en ΝΦΕΝ - 0 and 
Bar step is a ced. 

fereeach mid-side node that is on a curved element edge 
sS prepared with format (15,2F10.5) to read the global 
node number and the R and Z coordinates. Here, as in step 5, 
E cards may be put together in any arbitrary order. For 
τ. 6ος Ε1ρ.]16, there will be only one card to be punched 


Endet would read: 


{ων د‎ о E cT 
A AO LE 


Seer cm ° б ` edi qM s 


20500000807 0000060 0000066000560000000006666566500000202021C602í 


: lI 12 13 14 15 18 17 18 19 20 21 22 2223 25 26 22 26 29 3231 32 3.34 25 3632 3: 32 40 11 4243 44 A5 46 4; A8 49 50 51/52 53 54 5€ 5 57 56 53 GOEL ¢ 
τ... ης; 555“... е ы е өе I 4 q + 4 i | | 


NITE TITTIES CS 


Mero assure in all cases that the double precision capability 
of the program is not compromised by less precise floating 
pont inputs, it is recommended that all such inputs be made 
in D-FORMAT. The F-FORMAT instructions merely allocate 


ΠΡ card fields. 
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Step 8: 


lites properties of the different materials are specified 
Sims Step. For each material a card must be prepared 
ME ues the modulus of elasticity (E), coefficient of 
ο πα expansion (AL), Poisson's ratio (PfI), thermal con- 
Ига та ту (TK), density (DENS) and specific heat (SHT). 


INC read: 
Ἕν, DD ОТК, DENS, SHT 


ene format (2D12.4, 4F8.3). The first card will be for 
material number 1, the second card for material number 2, 


etc. 


ep 9: 

ο this step a single card must be punched, Starting in 
ean one, which reads the word BRITISH or METRIC in accor- 
dance with the system of units used. 

Brezeoımletes the input of geometric and material prop- 


ου data. 


peep TO: 

Por a transient temperature problem only (no stress cal- 
culations) leave a blank card for this step and proceed 
ΠΠ то step 16. 

oT Stress problem or thermal Stress problem, in this 
Seepewe specify the type of problem and the structural 
boundary conditions. 

The following quantities, when pertinent, are to be 


specified. 
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ØMEGA Ш r edn тор тоо, about thE Z2 axis, in 
ПОМОГ тое рет minute 


PRES The external pressure applied to a boundary 
SC ent 
TEXT The number of known temperature vectors for 


ο αυ ο ο ο Oi thelial stresses 1s desired 


LOAD ШЕКИ атон ον all additional κπονπ ἱσαά 
vector. If there is one, L@AD = 1, otherwise 
LOAD = 0 

NNLT Total humber Of modes fixed in the longitudinal 
direction 

NNRT tal nuüunmnber of nodes fixed in the radial 
direction 


ΕΞ An indication for the plane-end boundary condi- 
non Τι desired IPLANE, =. 1 


EN cof these items is not applicable, the correspond- 
wM τσι is left blank. 
Now we have to punch a single card for this step which 


eads: 
ÜMEGA, PRES, IEXT, LOAD, NNLT, NNRT, IPLANE 


ο] οτπας (ZF10.4, 515). 


авер 11: 

Е пете 1S no pressure loading (PRES=0), omrt this 
ПИ ror non-zero pressure, the total number of nodes on 
pressure loading boundary segment  (NPNT) and the global 
node numbers of these (pressure) nodes (NPN(I)) must be 
БШ ееп ете. Only one segment is permitted. 


We read in: 
NPNT, (NPN(I),I-1,NPNT) 


with the format (ШО 


14% 





iN example, Fig. 16, if there is pressure inside the 


Pedy, then the data card for this step will be: 


s i 4 b 11 14 j 


ro 
EA 


Ἕνα 
а 
ρω 

= 

t 
er 

4» 
= 
mo 
d 
ο) 
— 
— 
2-5 
со 
-— 
uo 
ns 
ce 

ο 
ws 
πο 
δω 
rv 
ta 
ρω 
de 
`x 
«rı 
“3 
D 
-- 
α; 
ρω 
^ eo 
^» 
2 
^A 
a 
to 
^ — 
La» 
ек5 
ta 
b > 
La 
Фф. 
c 
` 
C 
e 
ч 
~ 
=d 
oo 
Ca 
< 
P 
c 
im 
ES 
Lu 
~ 
ль 
[25 1] 
Lu 
de 
m 
ыл 
Ф» 
ὧν 
de 
-: 
í 
> 
~ 
GE 
ore 
> 
ase 
€ 
un 
m 
t^ 
< 
L 
Lu 
un 
on 
n 
< 
ar 
-- 
in 
xs 
Lu 
чэ 
m 
e 
o» 
БЕЖ 
o» 


*" - -^ » + © .» © © 


E tep 12: 

EM re τς no additional load vector (L@AD = 0), omit 
pats Step. 

For LOAD = | read in the additional load vector with the 
imal (6F10.4). The components of the load vector are 
4 in the order: R component at node 1, Z component 
erde), R component at node 2, etc. I.e., READ 
(F(I),I=1,NDF), where NDF (the number of degrees of freedom) 


is equal to two times the total number of nodes (2*NNT). 


Step 13: 

imeem, estress problem, or thermal stress problem, there 
must be at least one node constrained against longitudinal 
motıon, 1.€., NNLT > 0, so we must specify here which nodes 
are to be constrained against longitudinal motion. These 
global node numbers NNL(I) are read in with the format 


О) і.е 
READ: (NN(I),I=1,NNLT) 


In our example, if global nodes 1 and 21 are fixed in the Z 


direction, we have a card that reads: 
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x 4 
Ф. a - t » > ۹ а 
- 9^ 
er — 3 
a a = = T p © - 


Settee Este urbe ares teva u νο 


000000000000000000000000000000005000000000002690200000089200501 
|o 4551895 


| - 4647 5 C 56515659 eC 61 E 
3 26 27? 28 23 30 31.32 23 34 35 36 37 33 39 40 41 42 43444, 46 47 48 49 50 51 57 53 54 55 í 
E Е Р E Б H A р a : а 2 n A : 2 r а ; 9 τ... 399999 4 9 q EP Ἢ 1035391 7 5 10107151 


ч ° υ "νε 


ENSE 14: 

Е there are no nodes fixed in the radial direction 
monet = 0), omit this step. 

For NNRT > 0 we must read in the global numbers (NNR(I)) 
Of those nodes which are to be fixed against radial motion. 


E CEormat 15 (1015), i.e., 
READ: (NNR(I),I=1,NNRT) 


wmode and 21 of our example are also fixed in the R 
cron, then the card for this step would read exactly 


Bene one in Step 15. 


ECB 15: 
¡there is no plane-end boundary condition (IPLANE = 0), 


omit this step. 

ГОТ the case of end-planes-remain-plane boundary condi- 
tion, i.e., (IPLANE = 1), we have to specify the total number 
s ules Of the right-hand end (NNRE) and the global node 
numbers of this end (ЧҮК(1)). 


We read in: 
ШЕК СТ) ТЕМИКЕ) 


I Il tlle format (1015). 
rhe Case of Fig. 16, 1£ 1t 1s desired to have end- 
planes-remain-plane boundary condition then the data for 


Besten will be: 
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x tax s> » as j » 

Annan; 

0000000000000 00000000000000000000000080000000C00C20000C080CD6l 

1234567989101 12131416 16 12 18 18 2021 22 23 24 28 26 21 28 29 20 27 32 33 34 353 3€ 91 32 22 0 € 42 4344 48 46 4/ 49495051 52 93 54 35 55 57 $6 35 1U οἱ ὁ 

a a aaa q q q $ q q q q qa e Ú> q $ q 4 q q q <ç л у A + 454 4.4 ο ο ο τσ το 34 * 1 1 1* 3 3 
Sep 16: 


MENS step we begin the specification of thermal 
ΙΙ conditions. Additional details are prescribed in 
meeps 19 through 24. 

Thermal boundary conditions are imposed along discrete 
segments of the boundary. Each such segment must begin and 
end at a corner node of a boundary element. We consider 
Separately the imposition of the various kinds of thermal 
boundary conditions. 

fay Convection 

Simce the fluid temperature for convection is deter- 
Mined from a prescribed temperature-time history at entry 
"τοπ πα a specified flow velocity (see Appendix D, 
Perron 5), it is necessary that the terms "inside," i.e., 
adjacent to the symmetry axis, and "outside" have their 
wna meanings. Thus on the inside surface the convec- 
pron poundary condition may be applied to a single (con- 
menos) seement. A similar prescription may be employed 
ene Outside portion. The entry section used for flow 
ЕИО опо is at the upstream end of corresponding segment. 

(b) Constant temperature 

ПОСТ Ру constant temperature on a portion of the 

boundary, the designator "inside" or "outside" may be used. 


Such a constant temperature portion may consist of several 
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ШОО С с segments. If two different constant temperature 
ES S ire prescribed, one may be designated "inside" and 
mae Other "outside." 


fey Combinations of convection and 
"omscdgmtortemperature 


The convection and constant temperature conditions 
K used together, but the portion designated "inside" 
must have only one of these conditions prescribed and the 
mamemrestriction applies to the "outside." 

(d) Insulated 

i portions of the boundary not included in the 
S IS Specifically identified as "inside" and "outside" 
are considered insulated. 

The following items, when pertinent, are to be 


Piven as specified below: 


MEN TT The constant initial body temperature 
Q2=0. For insulated boundary condition outside 


Se 0s) For convection or constant temperature boundary 
condition outside 


Q3=0. For insulated boundary condition inside 


oe meer convection or constant temperature 
boundary condition inside 


TE It solving stress problem only 


INTE solving temperature problem only 


Q4=0. le olivine thermal stressi probiem 

КОО ЕТИ с ахат force in the Z direction. As usual, 
РИО TOT tension and MINUS] for compres- 
Sion. ° 


eo, we read: 


MIT 02, 03, O45 AXIALF 
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menthe format (5F10.4). 

GO ο ου το, ]6 TÉ the initial solid temperature is 
60° and we have insulated outside and convection boundary 
Exon inside with 120 kg axial compressive force, then 


ШИЕ ата for this step will bo: 


κ. 


00 000000009000000 0000006250000000C6 t6 G066GCOCQQOCSCÓLl 
E 


19 11 12 1314 151617 18 19 20 21 22 23 24 25 25 27 28 29 30 31 32 32 34 3, 35 37 39 ч о ἐ] 48 49 5, 5; $2 53 54 52 56 57 38 59 o0 61 6 
a ooo o 0 0t 9g o 50 O0 099999 9 9 <q < 


са ае а ОТОС еее Уе з 1 "SESS 1: ται 1 373-4 * 1 


Dnus the blanks left in the columns 51 through 40 indicate 


Bade want to solve a thermal stress problem. 


Sep 1!: 


Nor action 15 taken in this step - we merely choose between 
mesedrne to step 18 or jumping το ους ο 

If there are no known temperature vectors for which 
mation Of thermal stresses ais desired, proceed to 
EN Cp 18. 

EDO ies) 1.G., when some known temperature vectors 
are to be entered for thermal stress analysis alone or 
combined with some other loadings, proceed directly to 


s 5726. 


peep. 18: 


If (04 <0), i.e., we are to solve only a stress problem, 


proceed directly to step 27. 
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Seen 19: 
ШО στο iS an insulated boundary condition outside 
ЕЕЕ 0); omit the following steps and start from step 22. 
КОО ые Case Of a convection or constant temperature 
πας condition outside (Q2»0.), we have to specify 
the outside Medi emaInsSterseoOctricClene (HICA (for constant 
temperature HTCØ = е Е eM cons τοι πας Initial 
uu temperature (TEMP0) which may be equal to the solid's 
initial ~ up rature, the total number of nodes in contact 
with fluid outside (NCFQT), and the number of temperature 


ramps for outside flow (NRAMPQ). We read: 
HTCÓ, TEMPO, NCFÓT, NRAMPÓ 


Si ie format (D16.4,F10.4,215). See example in step 22. 


ошер 70: 

Here we read in the global node numbers of the nodes 
in contact with outside fluid (NCFQ(I)). The sequence of 
these node numbers must be in the direction of flow velocity. 


Read: 
(NCF@(I),1=1,NCFQOT) 
πμ {πο format (1215) (see example in step 25). 


Beep 21: 


οι οσους ramp of outside fllow we read in the time when 
АШЫ ΙΤ; (aris BDRYO(1,1), the initial temperature of the 
ramp BDRYO(1,2), the final temperature of the ramp (specify 


Ch I It differs from the initial temperature of the next 
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meno) BORYY(1,3) and the velocity of the fluid for this 


ramp BDRY#(1,4), with the format (4F10.4) so we read: 
((BDRY@(1,J) ,J=1,4) ,1=1,NRAMPQ). 


The number of cards prepared at this step is equal to the 


number of ramps. See example in step 24. 


Step D. 


If there is an insulated boundary condition inside 
ΜΙ proceed directly to step 25. 

Pom the case of convection or constant temperature 
boundary condition inside, Q3>0. We have to specify the 
ende heat transier coefficient (HTCI) (for constant 
inside temperature HTCI#107°), the con cont inside initial 
mamas temperature (TEMPI), which may be equal to the solid's 
mitral temperature, total number of nodes in contact with 
inside (NCFIT) and the number of ramps of the inside 


Flow (NRAMPI). We read: 
Pcie. PEMPI, NCFIT, NRAMPT 


Ше format (D16.4, F10.4,215). 

Pomeour example assume the time variation of entry temper- 
Meter tor the inside flow to be as in Fig. 17. The heat 
Anne tericostticient is 175.0 with initial inside fluid 
Neme ature 70.0. The card for this step is: 

em E T 3 


Ad «3541 , A ھب‎ dr ο va есе, Рл. т. е d æ si «а= -" € Ф è 


0000000 00 000000000002 00000900080 0660660660: 
i 2 2 $4350 51 52 52 56 85 56 97 S6 o GU BIE 


` ` =. 9 s e ο ος . ®» . s a % 9 + . 4 q @ ° ò> 4. 4 q q 4 5 4 - : 54 í 9*4 € 4 «4 * a4 4 q > f ç + 3 | 4 4 + 4 s 1 4 1 


im mm 
U 
T eras. 
a 
<> 
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Temperature 






300 


ОО, 
165. 


80 


100. 180. Time 


Figure 17. Entry Temperature-Time Variation. 


Шет détails of this temperature-time history are speci- 


IM un step 24. 


ερ ἐδ: 


Here we read in the global numbers of those nodes which 
ExNcontact with inside fluid, (NCEI(I)) with the format 
“τ. lt is important to read these hode numbers in the 
Nun etıon Of flow velocity. 

wommexample, if the inside flow of Pig. 16 is from left 


bt then for this step the card reads: 


0000000000000000000000006000000602000000000606660600000C00608! 
123456785 18 16 17 8 19 20 21 22 23 24 25 26 27 28 23 30 31 32 23 34 3S 35 37 14 33 10 41 42 43 CL 4S 43 4] 18 49 5051 5253 24 55 56 57 58 5950618 


€ < o @ @ « е oe @ @ @ @ @ ч ө «е «4 4 « GG q q æ 4 


ρω 
i 
δω 
coo 
ra 


VEAS «el. c ο ο «wc ο cw α еба AS 





peep 24: 


NIS Step data for inside flow are given. For each 
ie read: the time when the ramp starts BDRYI(I,1), 
ο. tial temperature of ramp I BDRYI(1,2), the final 
Ex usture of ramp I (specify only if. it differs from the 
ES! temperature of the next ramp)  BDRYI(I,3), and the 


EN ty or the fluid for this ramp BDRYI(I,4). We read in: 
ООШ jJ ЛЕ EJ КАМРИ) 


СЕЛИ отглас (4F10.4). 
umm ample of Fig. 17, let the flow velocity for the 
τ ΠΡ be 15 and that for the second and third segment 


bev20. Three cards are needed: 


lc BU Bn. n D ere ον 
ад 
ed 
3 CARD 
0000 00 0000 δῦ 00000 00 000062 09 6006006C05^^20200 
12345538 5101)]12131239 16 7 18 1926 0122 23 23 25 23 21 28 2:991 3« 3. .34 35 36 37 38 33 49 4] 42 5, 2345 4. 07 13 1050 5: 52" 
ino.n2 E ΠΠ ciu, np ΕΠ. ΤΩ 
nd 
me CARD 
000 00 00000000 000000. 00 8002 00 LU EEE te 
a E12 33 14 38 16 17 AU HI TE ¿1 22 2y 5 25 25 2 20 29 10 31 32 2 31 33 39 40 4! 47 23 41 4« tc € ip G51 £2 5s $4 $3 XE ο) 18 
nau nn, ng 15. Dn 
ST 
EC AI 
000 09 00000 00 0000000000000000520 0000000100 0000009 10000000 
1620374 S В 7 5 $S 1011121312 ΗΛ} 22 23 24 25 25 27 2229 39 3] a n pr $2 53 54 3S5 56 57 $8 5 59 δῇ 61 62 63€ 
Sue ren +q < EA 4 Αα ο ο ουσ αι л л 1} 1’ 
I p 25: 


lEEhiso step we must specify the following items: 


DTI INE time аСт ео (stopie) of trapezoidal 
integration 


TIMEI Memes hen ties brea lculated temperature 
ПОО ог 1S to Бе stored for thermal stress 
evaluation 
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EVERY The interval of time between the successive 
pcempcrsturesvectors whach are to be stored 


IVEC Total number of temperature vectors to be 
stored for thermal stress evaluation 


BNTP Ше aderon for printing -the calculated 
temperature vectors 


I£ INTP = 0 there will be no temperature prints 


πι TP 1 the program will prine the tem- 
Pr rature after every step of time integration 


If INTP = 5 the program will print the tem- 
πα ште atrter-every > steps or time Inte- 
νο εἴς. 


ENS step we prepare a single card that reads: 
INEST TMEIL EVERY; IVEC, ΤΝΤΡ 


with format (3F10.4,215). 

menesexample, if the step size of trapezoidal time inte- 
pon 1s 2 sec. and we desire to evaluate thermal stresses 
for 10 different temperature vectors, each 40 seconds apart, 
Ese rting with the first thermal stress calculation at 


Bu xcgonal to 60, then the card for this step is: 


r 


1 


£o veis 


-- 


e.nn ei, Di an. nn е 


000000 000000 00 0590000 00 0000000 0000000000 лате св 000000 ог 
23456783 
tere 3 1 11 


10 11 12 13 14 15 16 17 16 19 20 21 22 22 24 2$ 26 21 20 29 20 31 32 2? 34 38 36 37 38 38 40 41 42 43 01 25 1. 47 48 49 50 51 52 53 54 55 55 57 “5 4 606: 6 
1 το „ыл оа а 


11141911111 1 


1 η... τ, 4 1.9 1 2 1 3*2 оа за ο а з * 4. 43/079 ο ο ο $ q 


ο. 5n columns 34 and 35 indicates that the temperatures 
e printed every 15 steps of time integration. The 
"Ue of integration for this example is 420 seconds, 


Since 


ЮЕ = EVERY * (1VEC=1) = 420. 
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(op 26: 


Pi there are no known temperature vectors (IEXT = 0) 
hich thermal stress evaluation is desired, omit this 
Есер. 

Bor IEXT>0 we must specify the nodal temperatures 
Mie are to be stored for thermal stress evaluation. We 


read: 
EUSTOR (I,J) , J-1,NNT), I=1,IEXT) 


ОШ Т ре format (6F10.4). 

This means that we enter all the components of the first 
nodal temperature vector, followed by the components of 
the second nodal temperature vector and so on until IEXT 


such vectors have been entered. 


ορ 2T: 


The data cards for this problem are completed now. If 
no more problems are to be solved for the same geometry, 
nt 1 5 step. 

If another problem is to be solved for the same geometry 
and spatial discretization, increment by 1 the NPROB in step 


4 and start the input of the new problem from step 10. 


ошер ZO: 
I3 53s no access to the IBM 360 computer of the 


Naval Postgraduate School, omit the following steps and 
dene Tron step 32. 
For the convenience of the user, the author has put 


a so-called CHECK program and the program presented in 
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Appendix B in the computer system at N.P.S. (Naval Post- 
mm uate School). It is advised, however, that the user 
check his input data deck with the CHECK program before 
attempting to solve the problem. 

Bor Using the CHECK program prepare the following 


@eemtrol cards. 


//XXXX0000 JÓB (0000,0000FT,XX00) ,'NAME' , TIME-1 

//J@BLIB DD DSN=F0609.BAKH,DISP=SHR,UNIT=2314 , V@L=SER=DUFFY 
//GÓ EXEC PGM=CHECK,REGION=100K 

//FTOGF001 DD SYSQUT=A,DCB=(RECFM=FBA,LRECL=133,BLKSIZE=3325), 
/ / UNIT-SYSQUT | 


ZEXNUSFOOl] DD * 


Ku the first card is the regular FORTRAN job card used 
SAS institution. 
repare your deck as Fig. 18 and it may be read in from 


ШЕПОТ сата reader. 


|/* 


Data deck 


| Control cards oft step 28 
| Job card 


Bohrer ld. SEL UD о Usin: HEC Program. 


J. 





ντι Output of the CHECK program has any error messages 
or has not been run, then there are some mistakes in the 
Eum deck or control cards. 

If there are no error messages on the output of CHECK 
Em, then study carefully the output and make sure it 
is the same problem that has been intended to be solved. 
EX rhe correctness of the input data has been insured, 


“εδρα to the next step. 


осер 29: 


melo Stands for Axisymmetric Transient Thermal Stress 
and this is the name given to the program presented in 
Appendix B when stored in the computer. For using that we 


Mp repare the following control cards. 


//XXXX0000 JOB (0000,000FT,XX00) 'NAME' , TIME-10 

//JOBLIB DD  DSN-F0609.BAKH,DISP-SHR,UNIT-2314 , VÜL- SER- DUFT Y 
//GÓ EXEC PGM=AXITTS,REGION=425K 

//FTOGF001 DD SYSQUT=A,DCB= (RECFM=FBA,LRECL=133,BLKSIZE=3325), 
// UNIT=SYS@UT,SPACE= (CYL, (6,1)) 

//FTOSF001 рр * 


Mere the irst card is the regular FORTRAN job card. It is 
ОЕ edito ask for 10 minutes time, 1.e., TIME=10, since 
ο vould not affect the priority of the job within the 
class K Jobs. 

mirc the deck as in Fig. 19 (it may be read into the 


ООШ Г ет from the hot card reader) and submit a so-called 





τες request card, available in the computer center, to 


E operator on duty. 


/* 
Data deck 


Control cards of step 29 


Regular FØRTRAN job card (TIME=10) 


eure 19. Deck Set-up for Using AXITTS Program. 


een 50: 
¡the user wishes to have a listing of the program he 
i pare the following control cards as in Fig. 20 since 


Memo ram is also listed on the data cell for this purpose. 


Regular FORTRAN job card 
ЛИТИЙ EXEC PGM-IEBPTCH 
//SYSPRINT DD DUMMY 
//SYSUTI DD DSN=F0609.BAKH,DISP=LD,UNIT=2321 , 
7 VÓL-SER-CEL0032,DCB- (RECFM-FB,LRECL-80,BLKSIZE-2000) 
//SYSUT2 DD SYS@UT=A,SPACE=(CYL,6) 
A IN DD * 
PRINT TYPØRG=PØ,MAXFLDS=1 ,MAXNAME=1 
MEMBER NAME=AXITTS 
RECORD FIELD=(80) 
ps 


Г ште 20. Deck Set-up for Obtaining a Listing 
о the Programs. 1 is. 
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ΜΝ t up deck of Fig. 20 may be read in from the hot card 


сет to pet a listing of the program. 


шер 31: 


If the user wishes to obtain a deck of the program AXITTS 
πα prepare a deck as in Fig. 21 and read it in from the 
card reader. He must notify the operator on duty to be 


prepared for punching almost two boxes of IBM cards. 


Regular FORTRAN job card 
77/PUNCH EXEC PGM-IEBPTCH 
//SYSPRINT DD DUMMY 
//SYSUT1 DD DSN-F-0909.BAKH,DISP-QLD,UNIT-22321, 
7 VØL=SER=CEL003,DCB=(RECFM=FB, LRECL=80,BLKSIZE=2000) 
КОО UT? DD SYSØUT=B,SPACE=(CYL,6) 
SS IN DD * 
PUNCH TYPQ@RG=P@,MAXFLDS=1 , МАХМАМЕ=1 
MEMBER NAME=AXITTS 


RECORD FIELD=(80) 


/* 
Γον Όντο ο]. Deck Set up говото a Punched 
Deck of the Propram ΑκΙΤΤς 
жор 52: 


Isis ns the program presented in Appendix B, it the 
mrn Goce: not have access to the computer facility at N.P.S., 


he must punch a copy of the program. (Good luck!) 





mI τπο storage location 425K bytes are required and 
ου the Output may be longer than what usually is allowed, 
Eee additional πη. are recommended. The execution 
Eu required depends on the size of the problem and the 
un r of time integration steps; however, 10 minutes 
Eurer time would be sufficient for a problem of 140 nodes 
ОШОО steps of time integration, with 20 temperature 


EL rs stored for stress calculation. 
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APPENDIX “D 


PROGRAMMING 


M Computer program presented in Appéendix< B is formed 
from one main program and fifteen different subroutines. 
The communication between the main program and the sub- 
momemnes 1S handled through the common blocks. In order to 
minimize storage requirements, most of the storage location 
m System stiffness matrix is used for temperature 
@eweuiations. This is accomplished by use of EQUIVALENCE 
EE ments. The total storage requirement is 425K bytes. 

M section the function of some parts of the pro- 
gram is discussed and the assumptions used are brought to 


umtention. 


1. MAIN PROGRAM 
The main program is simply calling different subroutines 
when they are needed. The subroutines which are called in 


main program are: 


INPUT, PR@B, CANDY, FLOW, F@RMV, TEMPER, PL@T, STIFF, FQORMF, 


CENTF, PRESS, DISPL, AND STRESS. 


BEE SUBROUTINE INPUT 

In this subroutine the data regarding the geometry of 
N Mocs, subdivisions into elements, and the material 
Beoperties are read in and printed e: Calculation ot the 


mid-side coordinates based on the straight line is done 
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ПОО This subroutine calls for subroutine PLOTRZ once 

to produce a graphical representation of the nodal points 
then the half-band-width of the system stiffness matrix is 
RENculated from the connectivity array. Also, depending 
upon the choice of the system of units, the necessary con- 
yersions in each system of units are accomplished and the 
Memercnce temperature (70°F or 20°C) based on the system of 


pants 15 selected here. 


See SUBROUTINE PROB 

e ach problem this Subroutine is called once by the 
main program to read in and print out the information 
G ring the nature of the problem. Based upon the given 
information the type of the problem is distinguished here 
Er the thermal problems the length of time integration 


calculated. 


4. SUBROUTINE CANDY 

Exbnontine CANDY (C and Y) evaluates the capacitance 
matrix C and admittance matrix o ο M aie banded Torm. 
Por each problem this subroutine is called once. The func- 
muUa low chart of the subroutine CANDY is given in 
uM" w5ince the only non-zero elements of the symmetric 
matrix ү" η. [1] τε the diagonal eam απα second super: 
uMENPenads. in order to conserve storage and reduce the num- 
NG. Arithmetic Operations, the non-zero elements are 


ο j jn three different vectors (DIAG, QOFFl, @FF2). 
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L=1,NET au uU 


T=1,NGP CE=0, YE=0 





Evaluation rot 
oN oN 
> οξ > on 


οι πας του ο- 
Jacobian natrix YT 


Con truction Oi 


aN 3Ν 
GLE oR 


Calculation OT 
element capacitance and 
admittance matrices 


Construction ο Ө απο. η 
οσα σον ο. ү* 
were 


С CX DTE ا‎ 
С TEE 
















Performing Cholesky 
decomposition on new C 


ἵ “κε 22. Functional Flow Chart of CANDY. 
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ame che system capacitance and admittance matrices C 
* 
EX are formed, then the matrices A and G (Eq. 29) are 


KS I Cd and replace C and Y respectively. 


5. SUBROUTINE FLOW 

It there is any convection or constant temperature 
MEEEmal boundary condition, this subroutine is called from 
the main program for each step of time integration in order 
Eu 91 пате the temperature of the fluid nodes. The calcu- 
Neon 1S based upon the constant fluid flux assumption 
SS Doth inside and outside flow). Consider a section of 
BENE eprular cylindrical pipe as Fig. 23 and focus attention 


element 1+1 оп the inner boundary. 


Ὦ 


ντ 
N 


| 
| 
| 
| 
| 
| 
| 
| 


— - ------- : ———  - —— dE» 


2 


SS EE 
CIN 


ΠΠ», Fluid Node Representation for Inside Flow, 


aid iS moving from left to right we may number the R 
ο... coordinates of the fluid nodes of the element 1+1 as 
ШП 25, then the volume of fluid pumped into the pipe 


at time t is given by 
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Qin = a V T, (50) 


Ка апа у are, respectively, the area and velocity at 
E xentrance section. Also the volume of the pipe up to 
BI Emrid-node of element i*1 is 
2, 
ο τν sf mR dz. (57) 


21 


where M. 1s the volume of the pipe between the entrance 
Eon and element i*1 on the inner boundary. 


Now we can write 
3 
R= у В. N. (58) 


where N; are the one dimensional shape functions (given in 
pendix A). 
If we substitute Eq. 58 into 57 and make the coordinate 


Mumsftormation, after integration we get 


E ie 2 NN S _ | 
КЫ 21) (SIRT*64R)*RZ*40R]R,-8RjR.-14R,R.). (59) 


Now by equating Eq. 59 to Eq. 56 at any time r we may deter- 
mine whether the front of the flow has already passed the 
ПЕСО J Odo Of the wetted side of the element i+1. 

ο nq ar argument can be Carried ους tor the end node 
of the element itl. In that case we also get an equation 
ОЕ Л О μα. 50 wrth different numerical coefficients. 

tor the case Of the outside flow it has also been 


πη. το have a constant fluid flux and a fictitious 
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entrance circular cross-sectional area (whose radius is the 
largest R plus unity) has been assumed. 

lecnumerical coefficients of R's in Eq. 59 are used 
in subroutine FLOW and it has been assumed that the temper- 
ature of a fluid πιεί τν: mot cenando passage 
through the active section. This is believed to be an 
acceptable approximation for representative values of flow 
n ty and active length. With this assumption, for. a 
S (ry time-temperature relation (inside and outside 
Ens subroutine evaluates the fluid nodal temperatures 


at any time. 


6. SUBR@UTINE F@RMV 

every step of time integration this subroutine is 
ΙΙ DY the main program to form the vector v of the 
meee hand side of the discretized finite element (Eq. 13) 
ir the given thermal boundary conditions. The program 


mgr us self explanatory. 


7.  SUBR®UTINE TEMPER 

itemnmodal temperatures are calculated in this subroutine 
an he trapezoidal time integrations. Irons' correction 
ip lied here after every 10 steps of time integration. 
The temperature vectors selected for the stress analysis are 
stored. The transient temperatures will be printed out at 


the desired interval of time. 
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8. SUBROUTINE STIFF 

о апу stress or thermal stress problem subroutine 
STIFF is called once by the main program to evaluate the 
EDuEIness matrix of the system. Once the stiffness matrix 
of the system is calculated then the desired structural 
boundary conditions are applied and, at the end, the Cholesky 
Meeomposition is performed. The functional flow chart of 


Мол пе STIFF is given in Fig. 24. 


9. FORME 

EUEsubroutine FOÜRMF the thermal load vectors are calcu- 
lated for as many as (IVEC) given temperature vectors. 
These thermal load vectors are the columns of a rectangular 
DUX E. The provision is made that no matter how many 
memeenmature vectors are given (always IVEC < 20) the IVEC 
number of columns of the F m filled with the thermal load 
τσι and the last column of F will contain the load 
vector corresponding to the unit end displacement for zero 
axial force when the plane-end boundary condition is applied. 


The flow chart of subroutine FQRMF is given in Fig. 25. 


ШОО SUBROUTINE CENTF 

If the system is rotating about the axis of revolution, 
this subroutine is called once by the main program to eval- 
uate the centrifugal load vector. This load vector is 
always placed in the first column after the thermal load 


EIUS IVEC*l position) in E. 
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aw, NET 


Form decomposed D matrix 


Element stiffness 
matrix = U 


I=1,NGP 





Evaluation of shape 
function and 
ƏN ƏN 
DNE DD 


Form Jacopian 
Form B matrix 


Evaluate element 
оспе Ета 


Add element stiffness 
matrix fos b waren 


stiffness in banded form 





Copy part Of Salle te ера 
planes-remoin pl n, p r. 


| 
Application of boundary 
condi trons 
Decompose system 
stiffness matrix 


RETURN 





Pure z4. Functional Ουτε οἱ 5ΙΤΕΣ. 


jo; 








Evaluate shape functions 


oN oN 
and 3r » πῃ 


Construct la oblan 
Form B Matrix 


Evaluate temper MUTE 
at the Gauss point 


Evaluate в! р ες 


Form element load vector 


Add to the system 
load vector 






l II VE G 





Form last Vee EO somes 
for unit атр сеш пг 


RETURN 


Paure. 25. Functional TIO ο ο ο FORME. 
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IM SUBROUTINE PRESS 

If there is any pressure acting on the system, the 
w cure load vector is calculated in this subroutine and 
vector 1S placed in the column next to centrifugal load 
Dor. 1f any, otherwise it will fill the position desig- 


ν.μ tor centrifugal load vector in E. 


me SUBROUTINE DISPL 

ln -cdrsplacemeut vectors are found in “this subroutine. 
Ace the stiffness matrix is already decomposed in banded 
Bonn, then oe displacement vectors one after another are 
obtained by back and forward substitution. The principle 
of superposition is applied here and finally the axial 
ime, if any, is corrected for the plane-end boundary 


eendıtıon. 


ome SUBROUTINE STRESS 
For every problem this subroutine is called by the main 


param once to evaluate the stresses at the points 
1 


/ 3 
КОК аса and printed for each displacement vector. For 





(n=+1, E= + Jof each element. The transient stresses are 
each problem the maximum mean stress and the maximum octa- 
I sesrins stress, together with the corresponding 


ames and locations, are printed. 


Meee LIMITATIONS 
Insel process of the development of the program dis- 
КОСЫ со far, 1t has been intended that al of the calcula- 


Ben and storage of data occur in the computer without 
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using any external devices such as disks or magnetic tapes 
so as to be able to solve any sizable problem in relatively 
Short time. 

For this reason the following limitations (Table VI) are 
Sees. Orth which give an overall size of 450K bytes to the 


program. 


TABLE UT 


MAXIMUM VALUES FOR PROGRAM PARAMETERS 


NBAND Half band-width of the system stiffness matrix 66 


NMAT Total number of different materials 5 
МЕТ Total number of elements 40 
NNT Total number of nodes 149 
NCFOT To AI Number. of nodes in contact with 57 


outside fluid 


NCFIT Ше number of nodes in contact with | a 
πρ πας fluid 


NNRE Totalesnumber of nodes of right hand end 5 


NNLT Total number of nodes constrained against any 37 
Keneitudınal motion 


NNRT Total number of nodes constrained against any 37 
Peedi motion 


ЕЕС Пав number of temperature COIS stored for 20 
thermal stresses 


NRAMP@ Total number of ramps for outside flow 15 
NRAMPI Total number of ramps for inside flow 15 
NPNT tal number of pressures mode: E 
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